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ABSTRACT 

Motivation: One of the challenging questions in modelling biological 
systems is to characterize the functional forms of the processes that 
control and orchestrate molecular and cellular phenotypes. Recently 
proposed methods for the analysis of metabolic pathways, for 
example, dynamic flux estimation, can only provide estimates of the 
underlying fluxes at discrete time points but fail to capture the com- 
plete temporal behaviour. To describe the dynamic variation of the 
fluxes, we additionally require the assumption of specific functional 
forms that can capture the temporal behaviour. However, it also 
remains unclear how to address the noise which might be present in 
experimentally measured metabolite concentrations. 
Results: Here we propose a novel approach to modelling metabolic 
fluxes: derivative processes that are based on multiple-output 
Gaussian processes (MGPs), which are a flexible non-parametric 
Bayesian modelling technique. The main advantages that follow 
from MGPs approach include the natural non-parametric representa- 
tion of the fluxes and ability to impute the missing data in between the 
measurements. Our derivative process approach allows us to model 
changes in metabolite derivative concentrations and to characterize 
the temporal behaviour of metabolic fluxes from time course data. 
Because the derivative of a Gaussian process is itself a Gaussian 
process, we can readily link metabolite concentrations to metabolic 
fluxes and vice versa. Here we discuss how this can be implemented 
in an MGP framework and illustrate its application to simple models, 
including nitrogen metabolism in Escherichia coli. 
Availability and implementation: R code is available from the authors 
upon request. 

Contact: j.norkunaite@imperial.ac.uk; m.stumpf@imperial.ac.uk 
Supplementary information: Supplementary data are available at 
Bioinformatics online. 
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1 INTRODUCTION 

It is generally impossible to simultaneously measure the abun- 
dance of all the molecular entities making up biological systems. 
In gene expression assays, for example, we typically measure 
messenger RNA expression, but not the activity of transcription 
factors and/or the occupancy of transcription-factor binding 
sites. Similarly, in metabolomic analyses (Chou and Voit, 2012; 
Voit, 2013), key metabolites can be measured using, e.g. mass 
spectrometry or nuclear magnetic resonance quantification, but 
it is rarely possible to comprehensively quantify the metabolites 
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even within a single pathway. Typically, more interesting than 
metabolite and enzyme abundance are the fluxes through bio- 
chemical reactions and metabolic networks (Orth et ai, 2010; 
Schuster et ai, 1999). Fluxes, v = (v\, . . . , v m ) T , correspond to 
the rates at which molecules, x = (x\, . . . , x n ) J \ are turned over 
by the m reactions; regulation of fluxes in light of changes in 
environmental and physiological conditions is also intimately 
linked to cellular physiology. 

Although the fluxes are of central concern, they are hard to 
measure directly. Estimates for intracellular fluxes can be ob- 
tained by tracking products from isotope-labeled ( 13 C and 15 N 
metabolic flux analysis) metabolites through the metabolic net- 
work (Blank and Ebert, 2012; Zamboni, 2011). However, such 
an approach is restricted to a metabolically steady-state analysis 
and is not appropriate for capturing dynamical flux variations. 
Instead, theoretical analysis has often progressed by assuming 
stationarity of the metabolic processes, which in turn allows 
for characterizing the sets of steady-state fluxes under a set of 
suitable assumptions (Klamt and Stelling, 2003; Schwartz and 
Kanehisa, 2006; Voit and Almeida, 2004). Flux-balance analysis 
is the most popular example of this strategy, but it becomes 
questionable once the steady-state assumption can no longer be 
upheld. Furthermore, as more data on enzyme abundance be- 
come available, we should attempt to include such information 
and the impact on metabolic processes (Colijn et al, 2009; Rossell 
et al, 2013). 

Here we provide a new framework that allows us to model 
metabolic fluxes and their dynamics, and which deals with the 
missing data problem in metabolic analysis in a flexible and con- 
sistent manner. Gaussian processes (GP) belong to the armoury 
of non-parametric Bayesian methods and have been widely used 
to describe dynamical processes (Kirk and Stumpf, 2009) and to 
infer hidden states, e.g. transcription-factor activities (Honkela 
et al, 2010). In applications to metabolic modelling, parametric 
approaches can offer potentially incorrect representations of the 
underlying fluxes (Voit, 2013). The strengths of GP models arise 
from their non-parametric nature, which enables us to put priors 
directly on a function rather than on the parameters of a para- 
metric function. With a multiple-output GPs (MGPs), single GP 
framework can be extended to handle many outputs, enabling us 
to learn the unknown relationships between metabolic species. In 
turn, MGPs can be used to infill the sparsely sampled data 
(Boyle and Frean, 2004). This means that by using MGPs, it is 
possible to impute the missing data in between the metabolic 
measurements more efficiently. 

Here we develop a more general framework that uses so-called 
derivative GPs (Solak et al, 2003), which allow us to link 
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metabolite abundance, x (or concentrations) and fluxes v. This in 
turn enables us to also treat time course data on metabolites and 
monitor the changes that occur in fluxes, e.g. over the course of 
physiological responses, such as to changes in the environment 
(Bryant et ai, 2013). 

2 METHODS 

2.1 GP regression 

Gaussian process regression (GPR) can be applied to recover an under- 
lying dynamical process from noisy observations. A GP defines a prior 
distribution over all possible functions, and to specify a GP, we need 
expressions for the mean and covariance function that describe the be- 
haviour of the system output over time (Hay kin and Moher, 2010). Below 
we review the standard GPR methodology. 

In a typical regression problem, we connect inputs x and outputs z via 
functions, z = /(x), where x = (x\, . . . , x n ) and z = (z\, . . . , z n ) are con- 
tinuous ^-dimensional real-valued vectors. The observed values of the 
dependent variable, z, can be related to the independent variables, f(x) 
through, 

where e is a noise term, which is here assumed to be independent and 
identically distributed according to a Gaussian distribution, e ~ A/*(0, cr 2 ). 
In GPR, we place a GP (Haykin and Moher, 2010; McKay, 1998) prior 
over the functions f(x), f ~ QV, meaning that at any finite number of 
input points X\, ...,x„ the values flxi) have a multivariate Gaussian 
distribution with zero mean and covariance function, K, 

\Kxi), ...,f(x n )] T ~Af(0,K(x,x% 

Different functional forms can be chosen for the covariance function 
(Rasmussen and Williams, 2006), either to simplify computations or to 
reflect constrains imposed by the data. A flexible and generic choice is to 
set the covariance function to 

K(x p , x g ) = o\ exp ^- ^ |x^ - x g , 

where 9 = (a 2 , represent a set of unknown hyper-parameters, and x^ 
and x q are inputs. Thus, y = (y\, . . .,y n ) T has a multivariate normal dis- 
tribution with zero mean and covariance matrix C{6) = K + a 2 I, with I the 
identity matrix. The unknown set of hyper-parameters, 6, can be estimated 
from the data by evaluating the following log-likelihood function, 

C{6) = -\\o% \C{6)\ - \y T C(ey l y —Xogliz, (1) 

using either a maximum likelihood approach or by sampling from the 
posterior distribution with Markov chain Monte Carlo methods (Neal, 
1997). 

For any finite number of input (test) points, x\, . . . , x*, we define the 
joint prior probability distribution 

With the GP prior, it is possible to evaluate the posterior distribution 
over the functions; the values of / evaluated at inputs (x*, . . . ,x*) and 
conditioned on the observations y are jointly distributed as (Rasmussen 
and Williams, 2006), 

[/!#), ...J(x* r )] T \y~Af(m p ,K p ), (2) 

where 

m p = K(x* p , x q )[K(x p , Xg) + cr 2 l]~ l y, 

and 

K p = K(x* p , xp - K(x* p , x q )[K(x p , Xg) + cr 2 I]K(x p , xp. 



Although Equation (2) defines an appropriate GP posterior, which 
allows us to make predictions about a single variable y, it remains unclear 
how to deal with several variables simultaneously: if outputs are corre- 
lated then the standard GPR framework may fail in providing an ad- 
equate description. 

2.2 Multiple-output GPs 

Boyle and Frean (2004) introduced MGPs, where a set of dependent 
GPs is constructed via multiple-input multiple- output linear filters. This 
perspective can capture the dependencies between several variables by 
solving a convolution integral and specifying a suitable covariance func- 
tion, which in turn includes the cross and auto correlations among related 
variables. Our construction of derivative processes below builds on MGPs. 

Dealing with linear filters is central to signal processing where such 
filters describe a physical systems that can generate an output signal in 
response to a given input signal (Haykin and Moher, 2010; Roberts, 
2008). Linear filters are characterized by their kernel function (an impulse 
response) h(t), and the output z{t) can be expressed via convolution 
integral, 

z(t) = h(t) <8> x(t) = j h(z)x(t - z)dz, 

— oc 

where the symbol '(g)' denotes the convolution operator. To transmit the 
signal that has the mathematical properties of a GP, the kernel function, 
h{t) must be absolutely integrable, i.e. 

j \h(t)\dt<oo, 

— oo 

Then if the input X(t) is specified to be a Gaussian white noise process, 
the output process, Z(t), will also be a GP. 

Specifying a stable linear time-invariant filter with M white noise pro- 
cesses as inputs, X\(t), . . . ,X M (t), K outputs, Z\(t), . . . , Z K {f) and 
M x K impulse responses results in a dependent GP model (Boyle and 
Frean, 2005). A multiple-input multiple-output filter can thus be defined 
as 

m °° r 

z k(t) = J2 / h m iT)X m {t-T)dx, 

— oo 

where h m k(t) are kernel functions and Z^{i) is the /rth output. As dis- 
cussed previously, the observed variables might differ from expected vari- 
ables owing to the measurement noise, and we thus consider 

Y k (t) = Z k (t) + W k (t), (3) 

where Wkif) is a Gaussian white noise process with variance o\. 

Multiple-input multiple- output filters are able to capture the relation- 
ships between several variables Yt(t); in the model, these kind of depen- 
dencies are build in via shared input noise sources that enable the 
specification of valid covariance functions. For the sake of simplicity, 
let the impulse response be a Gaussian kernel, h m k(t) = v 777 ^exp{— \{t— 
^mk) 2 A m k\. Then evaluating the convolution integral leads to the follow- 
ing covariance function, 

m ~ 

C i} {d) = J2 / h m ix)h mj {x + d)dx 

(4) 

m=l V mi ' m .) I J 

where S = A mj (A mi + A m j)~ x A m j and d = t a — % is the temporal separ- 
ation between two input points, (see Boyle and Frean (2004) appendix for 
derivation and generalization to multidimensions). Constructing 
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intermediate matrices Cy permits the definition of a positive definite sym- 
metric covariance matrix C between K variables, 



Cki 



Ckk + cr 2 K l 



[NxN] 



Here N = J2i=i * s total number of observations, and N t the number 
of observations of variable i. Having defined the covariance matrix, we 
can use the log-likelihood, which has the form (1) for the inference of the 
hyper-parameters 9 = {v m k, iJ, m t,Amk}- Again, following Bayesian frame- 
work, we can use the results from the GPR section to evaluate the joint 
predictive distribution (2) for all outputs. Alternatively, for a particular 
variable i, predictions can be made using the appropriate marginal 
distribution, which is Gaussian, with mean m/(Y) and variance var/(^), 
given by 

m f (0 = k T C l y, 
\ari(f) = K-k T C-\ 



(5) 



where 



:[(C i j(t , -t J - l )---C i j(t'-tj, Nj )]. 



2.3 Derivative processes 

For a GP that is derived through a linear filter, Y(t) = h(t) <g> X(t) + W(t), 
where X(t) is a white noise GP, h(t) is a kernel function and W(t) is an 
additive noise, it is easy to formulate the expression of a derivative process. 
Taking a derivative of Y with respect to t, it is possible to obtain a new 
process U that is also a GP (Boyle, 2007), 



dt 



Y(i): 



oo 

/Is' 



h{t-T)\X{x)dx = g{t)®X{t\ 



Thus, it is possible to construct the derivative process by convolving a 
white noise GP X(t) with a derivative kernel function g(t). This definition 
enables us to consider derivative processes and the corresponding original 
processes as a collection of dependent GPs. This is true because the de- 
rivative processes and the original processes are derived from exactly the 
same input, X(t). 

To construct a dependent model for several related variables 
Y = (Y\, . . . , Y K ) and their derivatives U = (Ui, . . . , Uk), it is necessary 
to define a suitable covariance structure, which in principal arises from 
the initial covariance function (4). For example, for a set of four depend- 
ent outputs (two original and two derivative processes), the following 
equations can be applied to compute the covariances (Girard, 2004; 
Kirk, 2011; Solak et al, 2003), 



Autocovariance function of derivative process U t 

d 2 



vvc, 



and Uj 



VVdjid) = cov[- 



IdYi 
\dt 



dYi 




t=t ' & 


J 


between 


two 


dYj 




t=ta dt 


t=tb) 



dt a dtjj 



d 2 
dt a dth 



C u (d); 



Cij(d); 



Covariance between original process Y t and corresponding derivative 
process U t 



VCii{d) = cov[ Y it 



dY f 
dt 



dtb 



Cu(d); 



Covariance between original process Y t and derivative process Uj 

d 



VCij(d) = cov[ Y t 



Let R denote a block matrix, 



dYj 
dt 



dt 



■Cyid). 



R 



H 

-22 / 



L = R r , 



Cn C n DC n DC 1: 
C21 C22 DC21 DC22 

which describes the correlations between observations Y = (Y\, Y2) and 
their 'function' values Z = {Z\ , Z2), and corresponding derivative vari- 
ables U = (U\, U2) evaluated at any finite number of test points 
t\, . . . , t r . In a similar manner, let H denote 

C21 



H 



DC U 
\DC21 



DC l2 
DC 22 



DC n 
DC21 
DDC n 
DDC l2 



DC n \ 
DC22 
DDCn 
DDC 2l J 



where the Cy matrices contain the correlations between functions Z\ and 
Z 2 evaluated at a finite set of test points t\, ...,t r ; DCy the correlations 
between functions Z = (Zi,Z 2 ) and derivative variables U = (U\, U2) 
evaluated at the same test points; and finally, DDCy consists of auto/ 
cross-correlations between derivative variables U\ and U 2 . The matrices 
R, L and H are building components of the overall covariance matrix K, 
which is symmetric and positive definite, 



K 



(C + o 2 \ R\ 



At a finite number of input points t\, . . . , t r , the matrix K allows us to 
place a joint prior over observations Y, functions Z and derivatives U, 

[Yi,Y 2 ,Zi,Z2,Ui,U 2 ] ~ A/"(0,K). 

Evaluating a GP posterior 

[Zi, Z 2 ,Ui,U 2 ]|[Yi, Y 2 ] ~ Af(m post ,K post ), (6) 

where 

m post = L[C + a 2 iy l R and K post = H — L[C + <7 2 l] _I Y, 

enables us to make joint predictions for the original and derivative pro- 
cesses simultaneously. Alternatively, if there is no need to sample from 
the posterior process, we can use marginal Gaussian distributions to 
make predictions for individual output. The marginal distributions for 
output i and its derivative process at any input point f, 



m Yi (t*) = k Yi [C + (T 2 lY Y, 
muXO = k z ,[C + or 2 l] -1 Y, 
var F .(0 = k - k Yi [C + cr 2 l] _1 k^., 



(7) 



™r Ui (n = r ] -k Zi [C + cr 2 l] k T z , 

where m Yi is the mean of the original process, the mean of the de- 
rivative process, var^ the variance of the original process and var^ the 
variance of the derivative process, and furthermore 



- h,\) 



Equations (6) and (7) can easily be extended to make predictions about 
any number of variables. 



K = C^Q)) + or 2 , 


t] = VVCu(0) 








/ VC iX {f 




Cj\(t* — t\,Ni) 






ky, = 










\Ci2(t*-t 2 ,N 2 )) 




\VC a {f 
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3 APPLICATIONS AND RESULTS 

To demonstrate the performance of derivative processes, we 
consider two simulation examples — a system of two oscillating 
signals and a simple model of linear metabolic pathway — before 
turning to a more complicated metabolic process and, finally, 
some real metabolic network data. The derivative processes can 
be used to address the flux estimation problem from time course 
data. Here GPs describe the dynamics of metabolites, and 
the corresponding derivative processes capture the functional 
forms of the associated fluxes. Below, all examples were imple- 
mented using the free statistical computing platform R www.r- 
project.org. 

3.1 Oscillating signals 

A simple oscillating signal can be expressed as 
z(t) = A sin(W + 0), where A is the amplitude, co = 2nf the an- 
gular frequency and 0 the phase angle. This is a particularly 
useful example because it is easy to evaluate the performance 
of derivative processes, as the derivative signals have a known 
analytic form. We consider a simple system that consists of two 
oscillating signals, z\(t) and z 2 (t), 

zi(0 = sin(20, => z\(t) = 2cos(2f), 

z 2 (t) = sm(2t + ^, z 2 (0, =2 cos (2* + ^), 

with t e [0,47r]. To model real experimental measurements, we 
add random noise to the simulated trajectories, Y\(t) = 
zi(0 + 6i, Y\(t) = z 2 (t) + e 2 , where e 2 ;~ Af(0, OA 2 ); we have 
observations of both signals at regular time intervals, 
Di = {t u , iu^f 10 and D 2 = {t 2J , Y 2 j}?^ = 10. To build a 
single model that captures the relationship between the two sig- 
nals, we apply the dependent GP framework (3) (K = 2) on a 
combined dataset D = {D\,D 2 }; each signal can be expressed as 
a superposition of three GPs — two of which are constructed via 
convolution between a noise source and a Gaussian kernel, and 
the third one is an additive noise. We set parameters A t of each 
Gaussian kernel to be exp(/}) and noise levels to 
o\ = exp(^i), cr 2 = exp(??2), leading to a set of hyper-parameters 
0 = (vi,fi, /jL\, /jL 2 , rji, 7] 2 ), i = 1, . . . , 4. To build the model the fol- 
lowing priors are chosen: vufi ~ (1, 2 2 ), rjj ~ J\f(— 2, 2 2 ) and 
\ij ~ Af(0.5, l 2 ), j = 1,2; the maximum a posteriori (MAP) 
estimate 6 is determined using a multistarting Nelder-Mead op- 
timization algorithm (Nelder and Mead, 1965). Dependent GP 
posteriors (6) allow us to make joint predictions about both sig- 
nals and their derivative processes at any finite number of input 
points, and the resulting posterior processes are summarized in 
Figure 1. From these posterior processes, it can be seen that the 
mean behaviour of our model agrees with trajectories of under- 
lying noiseless signals, and to make predictions about derivative 
processes, it is enough to consider only samples from the original 
sinusoidal trajectories. 

3.2 Linear pathway 

Next we consider a linear metabolic pathway with two regulatory 
signals (see Goel et al. (2008) Supplementary Material for de- 
tails), which is summarized in Figure 3a. Here the flow from x x 
to x 2 is negatively regulated by metabolite x 3 , and x 3 increases 



(a) sin(2t) (b) sin(2t + ^) 




02468 10 02468 10 



Time Time 

Fig. 1. Predictions with MGPs model for two oscillating signals, (a and b) 
Dashed lines represent true behaviour of noiseless sin(-) trajectories; dots 
correspond to the noisy observations for both signals (data); solid lines 
are the mean behaviour of the MGPs model (predictions with original 
GPs); light areas correspond to two standard deviations at each predic- 
tion point, (c and d) Dashed lines represent true behaviour of noiseless 
cos(-) trajectories; solid lines show the mean behaviour of the MGPs 
model (predictions with derivative processes); light areas correspond to 
two standard deviations at each prediction point 

the transformation of x 2 into x 3 . A set of ordinary differential 
equations (ODEs) can be used to describe the dynamics of these 
two metabolites, x 2 and x 3 (x\ is the constant external input), 

■^1 Vmax o 5 

X2 = K m (l+^) + xr X2 ' X3 ' (8) 

X 3 — ~^ 2 ^3 

To apply the derivative process approach, we simulate the 
ODE model with the following parameter values 
(V max ,K m ,Ki) = (18.6819,9.7821,0.5992) and initial conditions 
x 2 (0) = 1, x 3 (0) = 1. In this model, the concentration of X\ is 
assumed to be constant and equal to 2. The dataset consists of 
selected points from simulated trajectories with added Gaussian 
noise J\f(0, 0.05 2 ). Again we combine the 'noisy' measurements, 
and fit the dependent GP model to make predictions about the 
original trajectories and their derivatives. To obtain a functional 
expressions for fluxes vi and v 2 we need to estimate a dynamical 
variations of metabolic, x 2 ,x 3 , derivatives. The derivative pro- 
cesses provide the predictions for the left side of Equation (8) at 
any finite number of time points, whereas the original GPs de- 
scribe the solution on the same ODE (8). This enables us to link 
the metabolite measurements to metabolic fluxes. Figure 2 illus- 
trates the predictions with posterior processes, where solid blue 
lines correspond to the mean behaviour of the model, dashed 
lines to the original x 2 and x 3 trajectories and solid green lines 
to their derivatives. In addition, if we assume that we are able to 
measure flux V3 = x 3 ' 5 , we can obtain the functional expressions 
for fluxes vi and v 2 that are summarized in Figure 2c and d. The 
dark pink lines illustrate predicted fluxes from noisy metabolite 
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(a) 
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Time 
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(a) 



LINEAR PATHWAY 



(b) 



► x 2 x 3 

v 2 v 3 



BRANCHED PATHWAY 



(c) Flux v 1 



(d) 



Flux v 2 



10 15 20 25 

Time 



10 15 20 

Time 



Fig. 2. Predictions with MGPs model for linear metabolic pathway, 
(a and b) Dashed lines represent a simulated x 2 and x 3 trajectories 
from ODE model; dots correspond to the sparse noisy observations for 
x 2 and x 3 (data); solid blue/green lines are the mean behaviour of the 
MGPs model (blue, predictions with original GPs; green, predictions with 
derivative process); light areas correspond to two standard deviations at 
each prediction point, (c and d) Dark lines are predicted fluxes, light areas 
correspond to the confidence region, and dashed lines represent true be- 
haviour of noise-free fluxes Vi and v 2 (calculated from ODE system) 



measurements, dashed lines are real fluxes (calculated from 
ODEs (8)) and light pink area corresponds to the confidence 
region. 

3.3 Branched pathway 

We now turn to an example of metabolic pathway that was ori- 
ginally proposed by Voit (2013) (see Example of actual charac- 
terization)', Figure 3b illustrates a schematic representation of a 
branched pathway with two regulatory responses, where x 3 in- 
hibits the conversions of x\ into x 2 , and x 2 positively regulates 
reaction v 4 . The following ODE model describes the dynamics of 
the metabolites that are involved in this pathway, 



: 0.05 - 

: 1.1*? 



X2 



x 3 = \Ax° 2 - 6 , 



-0.75 



l.l.vV\v ; " • z.x.v'.'-v:-- 



- l.b# 



(9) 



where x\,x 2 ,x 3 denote the metabolites. For a given pathway 
(Fig. 3b), the change in metabolite concentration can be 
described by the differences between incoming and outgoing 
fluxes. For this reason, we are able to obtain the following ex- 
pressions for fluxes v\,v 2 , v 3 and v 4 , 



Xi = vi - v 2 

X 2 = V 2 - V3, 

x 3 = v 3 , 



V4, 



V 4 = X\ + v 2 , 
v 2 = x 2 +i 3 , 
v 3 = x 3 . 



(10) 



These expressions define a system of linear equations that is 
underdetermined, as we have more fluxes to estimate than 




(c) 



NITROGEN ASSIMILATION 




aKG - a ketoglutarate 
GLU - glutamate 
GLN - glutamine 



Fig. 3. Pathway information, (a) A simple linear metabolic pathway; red 
and green dashed lines correspond to the inhibition and activation sig- 
nals, (b) Illustrates a branched pathway with positive (green) and negative 
(red) regulatory signals, (c) Illustrates a metabolic pathway in E.coli, here 
Vi, i= 1...4 denote the fluxes; aKG, GLU and GLN correspond to the 
metabolites; TCA is a short notation for the citrate cycle in E.coli 



available equations, and it cannot be solved using standard 
Gaussian elimination techniques. For this reason, additional in- 
formation is required to uniquely determine fluxes vi and v 4 . In 
this example, we will focus only on estimation of fluxes v 2 and v 3 
from available data rather than try to address a uniqueness prob- 
lem of vi and v 4 . 

The above ODE model enables us to generate simulated time 
course data using the initial conditions x\(0) = 4, x 2 (0) = 1 and 
x 3 (0) = 2. Next, we apply the dependent GP framework (3) 
(K = 2) on the combined dataset T) = {D\,D 2 }, where 
D\ = U2,/,*2,/}£r 2 ° andZ> 2 = {t 3 j,x 3J }f2^ 20 contains the meas- 
urements of metabolites x 2 and x 3 with added random Gaussian 
noise «/V(0, 0.0 1 2 ) (we chose a low noise level so that predictions 
with derivative processes could be easily compared with the ori- 
ginal fluxes in the example in Voit (2013). For a set of model 
hyper-parameters 0 = (yufu Vi,r]2,t*'),i= 1, . . . , 4 we use the fol- 
lowing priors, Vi - (2, 2 2 ), f x - (-3, 2 2 ), rjj - N(-2, 2 2 ), 7=1,2 
and \x — A/"(0.5, l 2 ), and calculate the MAP estimate 6 as before. 
Figure 4 illustrates the predictions with posterior processes using 
Equation (7); (a and b) graphs summarize metabolite data. The 
dark blue lines correspond to the mean behaviour of the original 
GPs and agree well with simulated x 2 and x 3 dynamics; the green 
lines describe the derivatives of the same metabolites and can be 
understood as a slope estimates. In Figure 4c and d, dark pink 
lines illustrate the predicted metabolic fluxes v 2 and v 3 under 
consideration of pathway Figure 3b. From ODE model (9), we 
can calculate original fluxes over the time (in real situations this 
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Flux vs 




Fig. 4. Predictions with MGP model for a branched metabolic pathway, 
(a and b) Dashed lines represent simulated x 2 and x 3 trajectories from 
the ODE model; red dots correspond to the sparse observations for x 2 
and x 3 (data); solid lines are the mean behaviour of the MGPs model 
(blue, predictions with original GPs; green, predictions with derivative 
process); light areas correspond to two standard deviations at each pre- 
diction point, (c and d) Dark lines are predicted fluxes; dashed lines rep- 
resent true behaviour of fluxes v 2 and v 3 (calculated from the ODE 
system) 



would not be possible). Figure 4c and d shows a good agreement 
between predicted and original fluxes. 



3.4 Escherichia coli nitrogen assimilation 

Finally, we apply our technique to the experimental data from 
E. coli, where we have measurements of the abundance of several 
key metabolites involved in nitrogen assimilation. Nitrogen is 
one of the key chemical elements that acts as a nutrient for the 
cells; ammonium is a preferred source of nitrogen for E. coli 
growth (Schumacher et al, 2013; van Heeswijk et al, 2013). In 
E. coli, ammonium can be absorbed via two pathways: glutamate 
dehydrogenase (GDH) that operates during cell growth in am- 
monium-rich environments, and glutamine synthetase-glutamate 
synthase (GS-GOGAT) that operates during cell growth in low- 
ammonium conditions (van Heeswijk et al, 2013). Here, we are 
focussing on experimental conditions, where after a period of 
nitrogen starvation, the bacterial cultures are spiked with ammo- 
nium (Schumacher et al, 2013); Figure 5a shows experimentally 
obtained measurements for a-ketoglutarade (aKG), glutamate 
(GLU) and glutamine (GLN) metabolites over the time after am- 
monium spike; red dots correspond to a wild-type (WT) E. coli 
metabolic measurements, and in squares — isogenic glnG deletion 
(AglnG) measurements. Below we focus on the pathway 
summarized in Figure 3c, which includes both GDH and GS- 
GOGAT. For modelling purposes, we assume that fluxes v 3 and 
v 4 can be summarized by the overall flux v 3 that describes the 
flow from GLU to GLN, as there is not enough information to 
discriminate between them. From the pathway, we can construct 
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Fig. 5. Predictions with MGPs model for E. coli (WT and AglnG). (a) 
The symbols indicate experimentally measured concentrations of aKG, 
GLU and GLN metabolites (dots for WT, squares for AglnG). Solid lines 
correspond to the mean behaviour of dependent GPs model, (b) Predicted 
derivative behaviour for aKG, GLU and GLN metabolites, where solid 
lines correspond to the mean behaviour of dependent derivative pro- 
cesses, (c) Predicted fluxes v\, v 2 and v 3 for convenience, dotted line illus- 
trates horizontal 0-axis 



a system of linear equations that describe the dependencies be- 
tween fluxes and metabolites, 



aKG = V] — V2, 
GLU = V2 — V3, 
GLN=v 3 , 



vi = aKG + GLU + GLN, 
v 2 = GLU + GLN, 
v 3 = GLN. 



(11) 



We fit a dependent GP model (3) (K = 3) to the WT data and 
then to AglnG data (collected from a strain where glnG is 
absent). In the model, aKG is expressed as a sum of three GPs: 
the first GP describes aKG, the second expresses the relationship 
between aKG and GLU and the third one describes additive 
noise; GLN is modelled similarly. However, GLU is modelled 
as the sum of four GPs, where the first three describe GLU; 
the dependence between GLU and aKG; the dependence between 
GLU and GLN; and the fourth is an additive noise. Choosing 
kernel functions to be Gaussian hkif) = Vk exp{— ^t 2 Ak), we 
obtain the MAP estimate for all hyper-parameters (17 in total). 
The predictions with posterior process (7) are summarized in 
Figure 5, where solid blue lines describe predictions with depend- 
ent GP models for WT E. coli, and green lines for AglnG. Using 
the relationship (11), we can estimate fluxes v\, v 2 and v 3 
(Fig. 5c). 

To evaluate our predictions, we can compare flux v 3 and GS 
protein levels in WT and AglnG E. coli (see Supplementary Fig. 
SI). In E. coli, glnG encodes the transcription factor, NtrC (ni- 
trogen regulator) that controls GS expression levels, and in its 
active form, GS catalyses glutamine synthesis (van Heeswijk 
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et al, 2013). Experimentally, it was observed that in AglnG case 
protein, GS levels were significantly lower compared with the GS 
levels in WT E. coli (see Supplementary Fig. SIC and D). 
Because there is less enzyme available to catalyse the reaction 
in AglnG, the flux v 3 in the mutant will be noticeably reduced 
compared with the WT flux v 3 (see Supplementary Fig. SI A 
and B). 



4 DISCUSSION AND CONCLUSIONS 

Flux estimation has become central to many analyses into the 
metabolic processes and mechanisms. Typically, the estimates for 
a set of fluxes are obtained in a point-wise manner at discrete 
time points. It is clear that this fails to capture the temporal 
behaviour of the fluxes and additional consideration of paramet- 
ric models is compulsory to fully explain the fluxes; further, this 
approach is susceptible to noise that is present in experimentally 
measured metabolite data. 

Here we have addressed these problems and proposed a novel 
non-parametric Bayesian approach to modelling metabolic 
fluxes. This is based on MGPs that enable the construction of 
derivative processes. Because the derivative processes and ori- 
ginal processes share the same input source, we can complement 
the dependent GP model and make joint predictions about ori- 
ginal and derivative processes at any finite number of input 
points. Such derivative processes can be applied to characterize 
the temporal behaviour of metabolic fluxes from time course 
data — without having to make reference, e.g. transcriptomic 
data, to explain temporal variation — and here we have demon- 
strated the applicability on simple models and a real-world 
example. 

GPs, including our approach, propagate uncertainty in line 
with the assumed covariance structures. This can lead to 
large confidence intervals, especially if the dependencies among 
different observations are not considered explicitly. With increas- 
ing number of metabolic species within the pathway, the deriva- 
tive process approach might become computationally costly due 
to the inference of a large number of hyper-parameters and a 
matrix inversion step; however, this limitation potentially might 
be addressed by considering a sparse approximation for the full 
covariance matrix of all metabolic species (Alvarez and 
Lawrence, 2009). These can in principle deal with genome-level 
data. 
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